Sustavi nelinearnih jednadžbi

Problem. Nađimo rješenje ξ=(ξ1,ξ2,,ξn) sustava od n jednadžbi

f1(x)=0,f2(x)=0,fn(x)=0,

i n nepoznanica x=(x1,x2,,xn). Uz oznaku f=(f1,f2,,fn)T, ovaj sustav možemo zapisati kao

f(x)=0.

Opisat ćemo Newtonovu metodu i tri kvazi-Newtonove metode:

  1. Broydenovu metodu,

  2. Davidon-Fletcher-Powell metodu i

  3. Broyden-Fletcher-Goldfarb-Schano metodu.

Sve metode, uz zadanu početnu aproksimaciju x(0), generiraju niz točaka x(n) koji, uz određene uvjete, konvergira prema rješenju ξ.

Napomena. Opisi metoda i primjeri se nalaze u knjizi Numerička matematika, poglavlje 4.4.

24.7 μs

Newtonova metoda

Jacobijan ili Jacobijeva matrica funkcija f u točki x je matrica prvih parcijalnih derivacija

J(f,x)=(f1(x)x1f1(x)x2f1(x)xnf2(x)x1f2(x)x2f2(x)xnfn(x)x1fn(x)x2fn(x)xn).

Za zadanu početnu aproksimaciju x(0), računamo niz točaka

x(k+1)=x(k)s(k),k=0,1,2,,

gdje je s(k) rješenje sustava

J(f,x(k))s=f(x(k)).

Za računanje Jacobijana koristimo paket ForwardDiff.jl. Za crtanje funkcija koristimo paket Plots.jl.

12.0 μs
22.9 s
Newton (generic function with 2 methods)
78.6 μs

Primjer 1

(Dennis i Schnabel, 1996) Rješenja sustava

2(x+y)2+(xy)28=05x2+(y3)29=0

su točke T1=(1,1) i T2(1.18,1.59).

6.1 μs
f₁ (generic function with 1 method)
83.3 μs
500 ns

Nacrtajmo funkcije i konture kako bi mogli približno locirati nul-točke:

3.2 μs
7.8 s
406 ms
656 ms

Vidimo da su nul-točke približno T1=(1.2,1.5) i T2=(1,1). Štoviše, T2 je točno jednaka (1,1) (1 iteracija u trećem primjeru). Nadalje, metoda ne mora konvergirati (četvrti primjer).

3.9 μs
J₁ (generic function with 1 method)
25.5 μs
2×2 Array{Float64,2}:
 10.0  14.0
 10.0  -2.0
1.2 s
1234
5.1 s

Primjer 2

(Dennis i Schnabel, 1996) Rješenja sustava

x12x222=0ex11+x232=0

su točke T1=(1,1) i T2(0.71,1.22) .

6.3 μs
246 ms
17.4 ms

Primjer 3

(Dennis i Schnabel, 1996) Zadan je problem f(x)=0, gdje je

f(x)=(x1x22+x2ex31).

Točna rješenja su T1=(0,0,0) i T2=(0,1,0). Izračunat ćemo nul-točke s nekoliko početnih aproksimacija.

26.4 μs
J₃ (generic function with 1 method)
85.6 μs
1234
763 ms

Primjer 4

(Rosenbrock parabolic valley) Zadana je funkcija

f(x)=100(x2x1)2+(1x1)2.

Tražimo moguće ekstreme funkcije, odnosno želimo riješiti jednadžbu

gradf(x)=0.

30.9 μs
f₄ (generic function with 1 method)
30.6 μs
222 ms
1.1 s

Iz kontura vidimo da je primjer numerički zahtjevan, dok analitički lako vidimo da je jedina nul-točka x1=(1,1).

U ovom primjeru je vektorska funkcija f zadana kao gradijent skalarne funkcije pa Jacobijevu matricu računamo korištenjem funkcije FowardDiff.hessian() koja računa aproksimaciju matrice drugih parcijalnih derivacija polazne funkcije f.

6.9 μs
12
7
1.1 s

Primjer 5

Zadana je fukcija

f(x)=i=111(x3exp((tix1)2x2)yi)2, gdje su brojevi (ti,yi) zadani tablicom: i1234567891011ti012345678910yi0.001.01.04.12.21.25.21.12.04.01.001

Želimo riješiti jednadžbu

gradf(x)=0. Za razliku od prethodnih zadataka, gdje je kondicija κ(J)=O(10)

u zadacima (a), (b) i (c) i

κ(J)=O(1000)

u zadatku (d), u ovom zadatku je

κ(J)>O(106)

pa je metoda netočna i ne konvergira prema točnom rješenju x=(4.93,2.62,0.28).

7.5 μs
f₅ (generic function with 1 method)
2.9 ms
2.3 s
18.4 ms
42.0 μs
508 μs

Broydenova metoda

Za odabranu početnu aproksimaciju x0 i matricu B0, za k=0,1,2,, računamo redom:

Bksk=f(xk)(sustav)xk+1=xk+skyk=f(xk+1)f(xk)Bk+1=Bk+(ykBksk)skTsksk

Na ovaj način izbjegavamo računanje Jacobijeve matrice u svakom koraku. Možemo uzeti B0=J(x0), ali i neku drugu matricu.

23.4 μs
Broyden (generic function with 2 methods)
58.3 μs
775 ms
31.5 ms
85.7 μs
16.6 ms
40.4 ms
18.4 ms
34.9 μs

Davidon-Fletcher-Powell (DFP) metoda

DFP je optimizacijska metoda koja traži točke ekstrema funkcije F:RnR, u kojem slučaju je f(x)=gradF(x).

Za odabranu početnu aproksimaciju x0 i matricu H0, za k=0,1,2,, računamo redom:

sk=Hkf(xk)βk=arg minβF(xk+βsk)sk=βkskxk+1=xk+skyk=f(xk+1)f(xk)Hk+1=Hk+skskTykskHkykykTHkyk(Hkyk).

Za matricu H0 možemo uzeti jediničnu matricu, a za izvršavanje iteracije nije potrebno rješavati sustav linearnih jednadžbi, već se sva ažuriranja vrše s O(n2) operacija.

Jednodimenzionalnu minimizaciju po pravcu xk+βsk računamo tako što metodom bisekcije nađemo nul-točke usmjerene derivacije.

11.4 μs
Bisekcija (generic function with 2 methods)
54.5 μs
DFP (generic function with 2 methods)
102 μs

Primjer. Nađimo točku ekstrema funkcije

f(x,y)=(x+2y7)2+(2x+y5)2.

Funkcija ima minimum u točki (1,3).

7.3 μs
f₆ (generic function with 1 method)
56.4 μs
5
1.3 μs
g₆ (generic function with 1 method)
18.4 μs
658 ms
58.8 ms
50.4 ms

Broyden-Fletcher-Goldfarb-Schano (BFGS) metoda

BFGS je optimizacijska metoda koja uspješno traži točke ekstrema funkcije F:RnR, u kojem slučaju je f(x)=gradF(x).

Metoda je slična DFP metodi, s nešto boljim svojstvima konvergencije.

Neka je zadana funkcija F:RnR, čiji minimum tražimo, i neka je f(x)=gradF(x).

Za odabranu početnu aproksimaciju x0 i matricu H0, za k=0,1,2,, računamo redom:

sk=Hkf(xk)βk=arg minF(xk+βksk)sk=βkskxk+1=xk+skyk=f(xk+1)f(xk)Hk+1=(IskykTyksk)Hk(IykskTyksk)+skskTyksk.

Za matricu H0 možemo uzeti jediničnu matricu, a za izvršavanje iteracije nije potrebno rješavati sustav linearnih jednadžbi, već se sva ažuriranja vrše s O(n2) operacija.

Jednodimenzionalnu minimizaciju po pravcu xk+βsk računamo tako što metodom bisekcije tražimo nul-točke usmjerene derivacije.

11.1 μs
BFGS (generic function with 2 methods)
115 μs
282 ms
48.8 ms
54.6 ms

Julia paketi

Prethodni programi su jednostavne ilustrativne implementacije navedenih algoritama. Julia ima paket NLsolve.jl za rješavanje sustava nelineranih jednadžbi i paket Optim.jl za nelinearnu optimizaciju.

7.3 μs
21.4 s
f₁! (generic function with 1 method)
30.8 μs
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [-1.0, 0.0]
 * Zero: [-1.1834670032425283, 1.5868371427230779]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 6
552 ms
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.5, 1.1]
 * Zero: [1.0000000000000002, 1.0000000000000002]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 6
44.2 μs
146 ms
184 ms
10.0 s
 * Status: success

 * Candidate solution
    Final objective value:     5.375030e-17

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 5.13e-09 ≰ 0.0e+00
    |x - x'|/|x'|          = 5.13e-09 ≰ 0.0e+00
    |f(x) - f(x')|         = 9.67e-17 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.80e+00 ≰ 0.0e+00
    |g(x)|                 = 2.10e-11 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    35
    f(x) calls:    102
    ∇f(x) calls:   102
1.2 s
 * Status: success

 * Candidate solution
    Final objective value:     3.572548e-12

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 1.87e+05 ≰ 0.0e+00
    |x - x'|/|x'|          = 8.70e-02 ≰ 0.0e+00
    |f(x) - f(x')|         = 6.03e-16 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.69e-04 ≰ 0.0e+00
    |g(x)|                 = 9.69e-09 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    41
    f(x) calls:    137
    ∇f(x) calls:   137
87.0 ms
200 ns